#!/usr/bin/python
#coding=utf-8

"""
------------
Programma per la stima dei parametri che caratterizzano il laser
------------

Cerca i file nella cartella ./dati/*dat
Per ognuno di essi effettua il fit con la funzione

([0]*P/2*Err(-sqrt(2)*(x - [2])/[1]) + [3])*R   [V]

[0] = responsività del fotodiodo [mA/mW]
[1] = w(z) [mm]
[2] = posizione del centro del fascio [mm]
[3] = corrente di buio ib [mA]
Err(y) = Integrale da -inf a y di cost*exp(-y^2)
P = potenza del laser [mW]
R = resistenza di carico [kO]

***Sarebbe meglio non lasciare [0] e [3] come parametri liberi*** 

Ottiene w(z) per ogni z considerato
----------

I w(z) vengono interpolati da

[0]*(1 + ((x - [2])/[1])^2)^.5

[0] = w_0 [mm]
[1] = z_0 
[2] = punto di minima dispersione del fascio

"""

from __future__ import division, print_function
from ROOT import TF1, TCanvas, TGraphErrors, TMath
import glob
from math import sqrt

R = 22.8+3
P = .5e-3
file_w_z_name = 'z_w_zerr_werr.dat'
file_w_z = open(file_w_z_name, 'w')
file_r_z_name = 'z_r_zerr_rerr.dat'
file_r_z = open(file_r_z_name, 'w')
file_ib_z_name = 'z_ib_zerr_iberr.dat'
file_ib_z = open(file_ib_z_name, 'w')

#Per ogni file ottiene w(z) e lo salva in file_w_z
data_path = 'dati/'
for file_name in glob.glob(data_path + '*.dat'):
    z = float(file_name.split('.')[1])
    n = float(file_name.split('.')[2])
    z += 1e-1*n
    z_err = 1/sqrt(12)
    graph_v_x = TGraphErrors(file_name)
    formula_v_x = '([0]*{0}/2*TMath::Erfc(-sqrt(2)*(x - [2])/[1]) + [3])*{1}'.format(P, R)
    func_v_x = TF1('func_v_x', formula_v_x, 0, 40)
    func_v_x.SetParameters(0.5, 0.33, 19.0)
    graph_v_x.Fit(func_v_x, 'q', '', 0, 22.1)
    graph_v_x.SetTitle('V(x), z = {0} mm;Posizione del coltello [mm];Ddp [V]'.format(z))
    graph_v_x.SetMarkerStyle(7)
    canvas_v_x = TCanvas('canvas_v_x', 'canvas_v_x')
    graph_v_x.Draw('ap')
    #canvas_v_x.SaveAs(file_name + '.eps')
    prob = func_v_x.GetProb()
    r = func_v_x.GetParameter(0)
    r_err = func_v_x.GetParError(0)
    ib = func_v_x.GetParameter(3)
    ib_err = func_v_x.GetParError(3)
    w = func_v_x.GetParameter(1)
    w_err = func_v_x.GetParError(1)

    line = '{0}\t{1}\t{2}\t{3}\n'.format(z, ib, z_err, ib_err)
    file_ib_z.write(line)
    line = '{0}\t{1}\t{2}\t{3}\n'.format(z, r, z_err, r_err)
    file_r_z.write(line)
    line = '{0}\t{1}\t{2}\t{3}\n'.format(z, w, z_err, w_err)
    file_w_z.write(line)
    raw_input()

file_w_z.close()
file_r_z.close()
file_ib_z.close()

#Interpola w(z) e ricava w_0 e z_0
graph_w_z = TGraphErrors(file_w_z_name)
graph_w_z.SetTitle('W(z);z [mm];w [mm]')
graph_w_z.SetMarkerStyle(7)
formula_w_z = '[0]*(1 + ((x - [2])/[1])^2)^.5'
func_w_z = TF1('func_w_z', formula_w_z, 0, 100)
func_w_z.SetParameters(0.32, 10, -100)
#graph_w_z.Fit(func_w_z, 'q', '', 0, 100)
canvas_w_z = TCanvas('canvas_w_z', 'canvas_w_z')
graph_w_z.Draw('ap')
canvas_w_z.SaveAs('w_z.eps')
raw_input()

#Andamento della responsività
graph_r_z = TGraphErrors(file_r_z_name)
graph_r_z.SetTitle('Responsività;z [mm];r [A/W]')
graph_r_z.SetMarkerStyle(7)
graph_r_z.Draw('ap')
raw_input()

#Corrente di buio ib(z)
graph_ib_z = TGraphErrors(file_ib_z_name)
graph_ib_z.SetTitle('Corrente di buio;z [mm];I_{buio} [mA]')
graph_ib_z.SetMarkerStyle(7)
graph_ib_z.Draw('ap')
raw_input()
